/**
This file is part of Gel2DE.

    Gel2DE is free software: you can redistribute it and/or modify
    it under the terms of the GNU Lesser General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    Gel2DE is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU Lesser General Public License for more details.

    You should have received a copy of the GNU Lesser General Public License
    along with Gel2DE.  If not, see <http://www.gnu.org/licenses/>.

	Copyright (C) 2012, University of Bergen
*/

#include "GelImage.h"

GelImage::GelImage(TiXmlElement* aGelImageNode, GelImageTransformation* aGelImageTransformation, vtkFileOutputWindow *aLogFile)
{
	const char* aFileName;

	m_logFile = aLogFile;

	m_parameterMap = new std::map<const std::string, double>();

	bool dataReaderSetup = false;

	TiXmlNode *gelImageChildNode;
	for(gelImageChildNode = aGelImageNode->FirstChild(); gelImageChildNode ; gelImageChildNode = gelImageChildNode->NextSibling())
	{
		TiXmlElement *gelImageChildElement = (TiXmlElement*)gelImageChildNode;
		const char* name = gelImageChildElement->Value();
		std::string strName(name);
		if(strName == "FileName")
		{
			aFileName = gelImageChildElement->GetText();
			SetupDataReader(aFileName);
			dataReaderSetup = TRUE;
		}		
		else if(strName == "Parameter")
		{
			TiXmlNode *parameterChildNode;
			for(parameterChildNode = gelImageChildNode->FirstChild(); parameterChildNode; parameterChildNode = parameterChildNode->NextSibling())
			{
				TiXmlElement* parameterChildElement = (TiXmlElement*)parameterChildNode;
				const std::string parameterName = parameterChildElement->Value();

				double parT;
				int outcome = parameterChildElement->QueryDoubleAttribute("value", &parT);
				if(outcome == TIXML_SUCCESS)
				{
					m_parameterMap->insert(std::pair<const std::string, double>(parameterName, parT));
				}
			}
		}
	}

	SetupResampler();
	if(aGelImageTransformation)
	{
		m_transform->SetMatrix(aGelImageTransformation->getMatrix());
		m_transform->SetOffset(aGelImageTransformation->getOffset());
	}

	SetupIterator();

	CalculateIntensityStatistics();

	//Necessary in order to load correct transformation
	//without displaying images.
	m_resampler->Update();

	m_fixed = false;
	m_enabled = true;

	SetNormalizationMode(NO_NORM, 1.);
	SetBlurLevel(0.);

}

GelImage::~GelImage(void)
{
	delete m_parameterMap;
}

void GelImage::SetReleaseDataFlags(bool aReleaseData)
{
	if (aReleaseData)
	{
		m_stdReader->ReleaseDataFlagOn();
		m_datReader->ReleaseDataFlagOn();
		m_stdReader->GetOutput()->ReleaseData();
		m_datReader->GetOutput()->ReleaseData();
	}
	else
	{
		m_stdReader->ReleaseDataFlagOff();
		m_datReader->ReleaseDataFlagOff();

		m_stdReader->Update();
		ImageType::SpacingType spacing = m_stdReader->GetOutput()->GetSpacing();
		spacing[0] = 1;
		spacing[1] = 1;
		m_stdReader->GetOutput()->SetSpacing(spacing);

		m_resampler->Update();
	}
}

void GelImage::SetupDataReader(const char* aFileName)
{

	//Should not read data yet!
	m_datReader = Reader::New();
	m_datReader->SetFileName(aFileName);

	m_datReader->GetOutput()->ReleaseDataFlagOn();
	m_datReader->Update();

	FILE *fp = fopen(aFileName, "r");
	if(fp)
	{
		std::string imsg(aFileName);
		imsg += " loaded.";
		m_logFile->DisplayText(imsg.c_str());
		fclose(fp);
	}
	else
	{
		std::string imsg("ERROR: ");
		imsg += aFileName;
		imsg += " not found!";
		m_logFile->DisplayText(imsg.c_str());
	}

	std::string strFileName(aFileName);
	std::string strBaseFileName = strFileName.substr(0, strFileName.length() - 4);
	strBaseFileName.append(" st.png");
	const char* stdstr = strBaseFileName.c_str();
	SetupStdReader(stdstr);

}

void GelImage::SetupStdReader(const char* aStdFileName)
{
	//Should not read data yet!
	m_stdReader = Reader::New();
	m_stdReader->SetFileName(aStdFileName);
	m_stdReader->ReleaseDataFlagOn();

	FILE *fp = fopen(aStdFileName, "r");
	if(fp)
	{
		std::string imsg(aStdFileName);
		imsg += " loaded.";
		m_logFile->DisplayText(imsg.c_str());
		fclose(fp);
	}
	else
	{
		std::string imsg("ERROR: ");
		imsg += aStdFileName;
		imsg += " not found!";
		m_logFile->DisplayText(imsg.c_str());
	}
}

void GelImage::SetBlurLevel(double aSigma)
{
	if (aSigma == 0)
	{
		aSigma = 0.01;
	}
	m_gaussianFilterX->SetSigma(aSigma);
	m_gaussianFilterY->SetSigma(aSigma);
	m_gaussianFilterY->Modified();
	m_gaussianFilterX->Modified();
}

double GelImage::GetBlurLevel()
{
	double blurLevel = m_gaussianFilterX->GetSigma();
	if (blurLevel < 0.05) return 0.;
	else return blurLevel;
}
int GelImage::GetNormalizationMode()
{
	return m_normalizationMode;
}

void GelImage::SetNormalizationMode(int aNormalizationMode, double aGlobalScale)
{
	m_normalizationMode = aNormalizationMode;

	if (aNormalizationMode == NO_NORM)
	{
		m_shiftScaleFilter->SetShift(0.);
		m_shiftScaleFilter->SetScale(1.);
	}
	else if (aNormalizationMode == MEAN_NORM)
	{
		m_shiftScaleFilter->SetShift(0.);
		m_shiftScaleFilter->SetScale(aGlobalScale/GetAvgIntensity());
	}
	else if (aNormalizationMode == MEDIAN_NORM)
	{
		m_shiftScaleFilter->SetShift(0.);
		m_shiftScaleFilter->SetScale(aGlobalScale/GetMedian());
	}
	else if (aNormalizationMode == ZSCORE_NORM)
	{
		m_shiftScaleFilter->SetShift(-1*GetAvgIntensity());
		m_shiftScaleFilter->SetScale(sqrt(GetRMSIntensity()*GetRMSIntensity() - GetAvgIntensity()*GetAvgIntensity()));
		m_shiftScaleFilter->Modified();
	}
}

void GelImage::SetupResampler()
{
	m_resampler = ResampleImageType::New();

	m_shiftScaleFilter = ShiftScaleImageFilter::New();
	m_gaussianFilterX = RecursiveGaussianImageFilter::New();
	m_gaussianFilterY = RecursiveGaussianImageFilter::New();
	
	m_shiftScaleFilter->SetScale(1.);
	m_shiftScaleFilter->SetShift(0.);

	m_gaussianFilterX->SetSigma(0.1);
	m_gaussianFilterY->SetSigma(0.1);
	m_gaussianFilterX->SetDirection(0);
	m_gaussianFilterY->SetDirection(1);

	m_shiftScaleFilter->ReleaseDataFlagOn();
	m_gaussianFilterX->ReleaseDataFlagOn();
	m_gaussianFilterY->ReleaseDataFlagOn();

	m_resampler->SetInput(m_gaussianFilterY->GetOutput());
	m_gaussianFilterY->SetInput(m_gaussianFilterX->GetOutput());
	//m_gaussianFilterX->SetInput(m_datReader->GetOutput());
	m_gaussianFilterX->SetInput(m_shiftScaleFilter->GetOutput());
	m_shiftScaleFilter->SetInput(m_datReader->GetOutput());

	m_transform = GelImageTransformation::TransformType::New();

	GelImageTransformation::TransformType::MatrixType yFlipMatrix;
	yFlipMatrix[0][0] = 1;
	yFlipMatrix[0][1] = 0;
	yFlipMatrix[1][0] = 0;
	yFlipMatrix[1][1] = -1;

	m_transform->SetMatrix(yFlipMatrix);
	GelImageTransformation::TransformType::InputPointType p;
	//This we could do when needed.
	ImageType::SizeType imgsize = m_datReader->GetOutput()->GetLargestPossibleRegion().GetSize();
	ImageType::SpacingType spacing = m_datReader->GetOutput()->GetSpacing();
	spacing[0] = 1;
	spacing[1] = 1;
	m_datReader->GetOutput()->SetSpacing(spacing);
	p.SetElement(0, imgsize.GetElement(0) * spacing.GetElement(0)/2.);
	p.SetElement(1, imgsize.GetElement(1) * spacing.GetElement(1)/2.);
	m_transform->SetCenter(p);

	InterpolatorType::Pointer interpolator = InterpolatorType::New();

	m_resampler->SetInterpolator(interpolator);
	m_resampler->SetTransform(m_transform);

	m_resampler->SetSize( m_datReader->GetOutput()->GetLargestPossibleRegion().GetSize());
	m_resampler->SetOutputOrigin(  m_datReader->GetOutput()->GetOrigin() );
	m_resampler->SetOutputSpacing( m_datReader->GetOutput()->GetSpacing() );
	m_resampler->SetOutputDirection( m_datReader->GetOutput()->GetDirection() );
	m_resampler->SetDefaultPixelValue( 100 );

	m_resampler->Update();

}

void GelImage::SetupIterator()
{
	ImageType::RegionType region = m_resampler->GetOutput()->GetRequestedRegion();
	m_iterator = new ConstIteratorType(m_resampler->GetOutput(), region);

	m_iterator->GoToBegin();

}

void GelImage::SetRoi(ImageType::RegionType aRegion)
{
	delete m_iterator;
	m_iterator = new ConstIteratorType(m_resampler->GetOutput(), aRegion);
	m_iterator->GoToBegin();
}

void GelImage::UseData()
{
	m_stdReader->GetOutput()->ReleaseData();
	m_resampler->SetInput(m_gaussianFilterY->GetOutput());
//	m_resampler->SetInput(m_datReader->GetOutput());
	m_resampler->Modified();
}


void GelImage::UseStd()
{
	m_datReader->GetOutput()->ReleaseData();
	m_resampler->SetInput(m_stdReader->GetOutput());
	m_resampler->Modified();
}

void GelImage::Update()
{
	m_resampler->Update();
}

void GelImage::SetEnabled(bool aEnabled)
{
	m_enabled = aEnabled;
}

void GelImage::SetFixed(bool aFixed)
{
	m_fixed = aFixed;
}

bool GelImage::isEnabled()
{
	return m_enabled;
}

bool GelImage::isFixed()
{
	return m_fixed;
}

const char* GelImage::GetFileName()
{
	return m_datReader->GetFileName().c_str();
}

void GelImage::Rotate(double aAngle)
{
	if(!isFixed())
	{
		m_angle += aAngle;
		m_transform->Rotate2D(aAngle, FALSE);
	}
}

void GelImage::Scale(double aScale)
{
	if(!isFixed())
	{
		const double* scale = m_transform->GetScale();
		double newScale[2] = {0., 0.};
		newScale[0] += aScale + scale[0];
		newScale[1] += aScale + scale[1];
		m_transform->Scale(newScale);
	}
}


void GelImage::SetCenterOfRotation(double aCenter[2])
{
	GelImageTransformation::TransformType::InverseMatrixType invm = m_transform->GetInverseMatrix();
	GelImageTransformation::TransformType::OutputPointType point;
	point.SetElement(0, aCenter[0]);
	point.SetElement(1, aCenter[1]);
	itk::Point<double, 2> p = invm*point;

	m_transform->SetCenter(p);
}

void GelImage::Translate(double aVector[2])
{
	if(!isFixed())
	{
		GelImageTransformation::TransformType::OutputVectorType vector;
		vector.SetElement(0, aVector[0]);
		vector.SetElement(1, aVector[1]);
		m_transform->Translate(vector, TRUE);
	}
}

BOOL GelImage::HasParameterValue(std::string aParameterName)
{
	std::map<const std::string, double>::iterator pair = m_parameterMap->find(aParameterName);
	return pair != m_parameterMap->end();
}


double GelImage::GetParameterValue(std::string aParameterName)
{
	std::map<const std::string, double>::iterator pair = m_parameterMap->find(aParameterName);
	if(pair != m_parameterMap->end())
	{
		return pair->second;
	}
	else
	{
		return 0;
	}
}

GelImage::ImageType::Pointer GelImage::GetImage()
{
	return m_resampler->GetOutput();
}

void GelImage::MoveToNextPixel()
{
	++(*m_iterator);
}

GelImage::PixelType GelImage::GetCurrentPixelRaw()
{
	return m_iterator->Get();
}

double GelImage::GetAvgIntensity()
{
	return m_AvgIntensity;	
}

double GelImage::GetMedian()
{
	return m_Median;
}

double GelImage::GetRMSIntensity()
{
	return m_RMSIntensity;
}

double GelImage::GetSumIntensity()
{
	return m_SumIntensity;
}

void GelImage::CalculateIntensityStatistics()
{
	double sum = 0;
	double rms_s = 0;
	int N = 0;

	ImageType::RegionType region = m_datReader->GetOutput()->GetRequestedRegion();

	int npixels = region.GetNumberOfPixels();

	std::vector<PixelType>* pixelvector = new std::vector<PixelType>();
	pixelvector->reserve(npixels);
	MoveToStart();
	
	while(!IsAtEnd())
	{
		PixelType val = GetCurrentPixelRaw();
		pixelvector->push_back(val);

		double Axy = val;//GetCurrentPixelRaw();
		rms_s += Axy*Axy;
		sum += Axy;
		N++;
		MoveToNextPixel();
	}
	std::nth_element(pixelvector->begin(), pixelvector->begin() + (npixels/2), pixelvector->end());

	PixelType median = pixelvector->at(npixels/2);

	double rms_ms = rms_s/N;

	m_RMSIntensity = sqrt(rms_ms);
	m_SumIntensity = sum;
	m_AvgIntensity = sum/N;
	m_Median = median;
	std::cout << " Sum " << m_SumIntensity << std::endl;
	std::cout << " RMS " << m_RMSIntensity << std::endl;
	std::cout << " Avg " << m_AvgIntensity << std::endl;

	char text[1024];
	sprintf(text, " Mean: %.2f Median: %.2f RMS: %.2f ", m_AvgIntensity, (double)m_Median, m_RMSIntensity);
	m_logFile->DisplayText(text);


	delete pixelvector;

	MoveToStart();
}

void GelImage::MoveToStart()
{
	m_iterator->GoToBegin();
}

bool GelImage::IsAtEnd()
{
	return m_iterator->IsAtEnd();
}

std::set<const std::string> GelImage::GetParameterTypes()
{
	std::set<const std::string> types;
	std::map<const std::string, double>::iterator i = m_parameterMap->begin();
	for(i; i != m_parameterMap->end(); ++i)
	{
		types.insert(i->first);
	}
	return types;
}

TiXmlElement* GelImage::ToXML()
{
	TiXmlElement* element = new TiXmlElement("GelImage");

	//Filename
	TiXmlElement* fileE = new TiXmlElement("FileName");
	TiXmlText* fileT = new TiXmlText(GetFileName());
	fileE->LinkEndChild(fileT);
	element->LinkEndChild(fileE);

	//Registration transformation matrix
	GelImageTransformation::TransformType::MatrixType matrix = m_transform->GetMatrix();
	TiXmlElement* matrixE = new TiXmlElement("Matrix");
	matrixE->SetDoubleAttribute("e00", matrix[0][0]);
	matrixE->SetDoubleAttribute("e01", matrix[0][1]);
	matrixE->SetDoubleAttribute("e10", matrix[1][0]);
	matrixE->SetDoubleAttribute("e11", matrix[1][1]);
	element->LinkEndChild(matrixE);

	//Registration transformation offset
	GelImageTransformation::TransformType::OffsetType offset = m_transform->GetOffset();
	TiXmlElement* offsetE = new TiXmlElement("Offset");
	offsetE->SetDoubleAttribute("e0", offset[0]);
	offsetE->SetDoubleAttribute("e1", offset[1]);
	element->LinkEndChild(offsetE);

	return element;
}
